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Summary 



By reducing core body temperature, T c , induced hypothermia is a therapeutic tool to pre- 
vent brain damage resulting from physical trauma. However, all physiological systems 
begin to slow down due to hypothermia that in turn can result in increased risk of mortality. 
Therefore, quantification of the transition of T c to early hypothermia is of great clinical 
interest. Conceptually, T c may exhibit an either gradual or abrupt transition. Bent-cable re- 
gression is an appealing statistical tool to model such data due to the model's flexibility and 
greatly interpretable regression coefficients. It handles more flexibly models that tradition- 
ally have been handled by low-order polynomial models (for gradual transition) or piece- 
wise linear changepoint models (for abrupt change). We consider a rat model for humans 
to quantify the temporal trend of T c to primarily address the question: What is the critical 
time point associated with a breakdown in the compensatory mechanisms following the 
start of hypothermia therapy? To this end, we develop a Bayesian modelling framework for 
bent-cable regression of longitudinal data to simultaneously account for gradual and abrupt 
transitions. Our analysis reveals that: (a) about 39% of rats exhibit a gradual transition in 
T c ; (b) the critical time point is approximately the same regardless of transition type; (c) 
both transition types show a significant increase of T c followed by a significant decrease. 

Key words: Bayesian inference; bent-cable regression; brain damage; cardiac arrest; grad- 
ual and abrupt transitions; mixture model; transition point. 
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1. Introduction 

Longitudinal data arise in many areas of study, where measurements taken over time are 
nested within observational units drawn from some population of interest. In particular, 
data showing a trend that characterizes a change due to a system shock are commonly ob- 
served over time in biological, medical, health and environmental applications. An example 
is an experiment on 38 rats (Reynolds et al. , 2007, also see Section [2]) conducted with an 



objective to collect information about the state of hypothermia and resuscitation strategy 
immediately after a 60% hemorrhage. In practice, hypothermia results in an initial increase 
in core body temperature, T c , before a decrease takes place. However, critically low T c may 



result in a breakdown in the compensatory homeostatic mechanisms (Connett et al. , 1986 



Rincon & Mayer, 2006). Therefore, timely resuscitation from hypothermia is of great clin- 
ical interest, which requires (i) the identification of the critical threshold at which T c starts 
to decrease, and (ii) the estimation of the decrease rate in T c after the transition. 

Figure [T] shows six of the 38 temporal profiles of T c in grey. They are selected to 
reflect the range of shapes of the 38 profiles. Overall, a similar type of trend is exhibited 
by all profiles - roughly linear incoming and outgoing phases are observed at the ends 
of each profile, with a continuous transition between phases. Some rats exhibit a gradual 
transition, while others, an abrupt transition. That is, we have samples potentially coming 
from two populations, labelled, G (gradual) and A (abrupt), respectively, according to the 
type of transition for the underlying T c trend. An exception in the figure is Rat 4, which 
exhibits neither an obviously gradual nor abrupt transition, but rather a seemingly linearly 
decreasing trend. There are only four of such profiles in the dataset, not adequate for 
hypothesizing an additional population. Such an investigation could be possible with a 
sufficiently large dataset. 

[Figure 1 about here.] 
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Accounting for the possibility of two well-defined populations, A and G, given as few 
as 38 rats, we develop a statistical framework for modelling these data with particular inter- 
est to address (i) and (ii) mentioned above, among related questions concerning therapeutic 
hypothermia. Our modelling approach is a substantial generalization of a special change- 
point model, the bent cable ( |Chiu et alj|2006| ). It provides flexibility with which inference 
for the type of transition for each individual is data driven, rather than pre-assumed as a 
specific type. |Chiu et aL\ ( |2006| ) and |Chiu & Lockhart| ( |2010| ) developed the bent-cable 
regression methodology and inference asymptotics to analyze a single data profile show- 
ing roughly three phases: incoming and outgoing, both of which are linear, joined by a 
quadratic bend (Figure [2ja)). As an extremely sharp bend reduces the bent cable to a bro- 
ken stick (Figure |2fb)), the former encompasses the latter as a limiting case. Although the 
model is parsimonious and appealing due to its simple structure and great interpretability, 
the authors pointed out that the segmented nature of the model may lead to poor asymptotic 
approximation in many practical settings involving finite samples. 



[Figure 2 about here.] 



Khan et aL\(2009) showed an extension of the bent-cable regression for longitudinal 



data by explicitly hypothesizing that the sample came from Population G only (henceforth, 
we will refer to it as Model G), but virtually no methodological details were provided. 
Their emphasis was on an atmospheric phenomenon that took decades to develop a clear 
temporal trend. In contrast, in this paper, we investigate two distinct physiological trends 
that manifest themselves within minutes. For this, we develop our flexible methodology to 
account for both gradual and abrupt transitions, for which Model G is a simpler special case. 

The piecewise linear (broken stick) model has been heavily utilized to describe a con- 
tinuous trend exhibiting at least one abrupt change over time, (e.g., Bellera et al. , 2008 



Hall et al\ |2003| |Kiuchi et al.\ |1995[ ). However, abruptness of change for all individuals 
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necessarily imposed by the broken stick is unrealistic for the hypothermia experiment, as 
demonstrated by Figure [T] Other quantitative methods that ignore temporal correlation and 
treat each subject separately are also unsatisfactory. This is because temporal correlation 
is typically substantial ( Reynolds & Chiu[ 2010[ ). Moreover, based on the individual model 



fits, any inferential statement for the underlying population is at best ad hoc. The use of our 
flexible bent cable framework presented in this paper advances insight into the important 
aspect of the quantification of the T c trend over time by (a) relaxing universal abruptness 
or gradualness through mixture modelling of piecewise linear and bent cable, (b) allowing 
proper inference at the population level pooling all individuals in a mixed effects longitu- 
dinal framework, which also (c) mitigates practical difficulties with modelling, such as the 
need for ignoring outlier individuals and for arbitrary truncation of data (as were necessary 
for Reynolds & Chiu 2010, see Section [2]). Note the practical significance of (a): if all 



individuals are assumed to exhibit the same type of transition, then shrinkage towards the 
population may force an observed profile resembling a broken stick (e.g., Figure[l£) to take 
on a bent-cable fit, and vice versa. Such bias will be demonstrated in Sections [4. 1 1 and |6| 
That is, due to shrinkage, the broken stick being a special case of the bent cable (Figure [2]) 
does not necessarily prevent biased inference for an individual rat. Thus, property (a) is 
crucial for proper inference in the context of therapeutic hypothermia. 



Although flexible modelling approaches such as penalized spline regression (Ruppert 



et al. , [2003) can also handle abrupt and/or gradual changes, the added flexibility in the 



shape of the fitted model can come at a cost of interpretability. In contrast, we use a model 
that is simultaneously flexible, interpretable and parsimonious at the population level (see 
Section [3]), and therefore potentially valuable to many scientific contexts. 

In Section [2j we describe the rat study and outline some substantive research questions 
about hypothermia therapy. In Section [3} we present our modelling framework to account 
for either type of transition through a longitudinal mixture model extension of the single- 



profile bent-cable regression technique. Additionally, an autoregressive process (|Box et 



al. , 2008 ) of order p (AR(p)), p > 0, is considered to approximate the within-individual 



autocorrelation structure. This is all constructed under a Bayesian framework (Section [4]), 
so that the concern about unsatisfactory performance of bent-cable asymptotics is irrele- 
vant. We then apply our method to the aforementioned rat data (Sections [5]) to address 
the questions regarding hypothermia therapy. In Section |6j simulations (a) demonstrate the 
importance of hypothesizing both Populations A and G for the rat model, and (b) illustrate 
that the flexible methodology can perform well with respect to the population regression 
coefficients even for a misspecified within-individual correlation structure, a fact which is 
taken into account in analyzing the rat data. We summarize our findings in Section [7J 

2. Data and Research Questions 

Neuronal damage is a common outcome for the survivors of cardiac arrest. Cardiac arrest 
generally leads to a decrease in the level of oxygen, a condition called anoxia which our brain 
can tolerate for up to 2 to 4 minutes ( |Krause et al.\ |1986| ); irreversible brain damage be- 



gins to occur thereafter. In fact, anoxic brain injury is the outcome of a complex process - 
ischemia and the subsequent reperfusion together cause enormous biochemical, structural 
and functional insults that lead to progressive cell destruction, multiorgan dysfunction and 
neural apoptosis ( |Negovsky}|1988p . Hypothermia can protect the brain and heart by atten- 



uating or ameliorating the deleterious temperature- sensitive mechanisms of that process. 

Effects of hypothermia on metabolism include a decrease in cerebral blood flow and 
brain volume, reduction of metabolism, diminution of intracranial pressure, and inhibition 
of glutamate release and other pathophysiological mechanisms (Rincon & Mayer} 2006). 



It protects tissue from ischemic damage through this process (Gordon, 2001J ). In contrast, 



when the body becomes very cold, all physiological systems begin to slow down, even- 
tually to the point that threatens survival. Treatment priorities in such situations include 
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prevention of further cooling and resuscitation (MartynJ 198 Therefore, the main re 



search interest lies in quantification of the transition of T c to early hypothermia. 

Motivated by the above, [Reynolds gf al] ( |2007| ) conducted the aforementioned rat exper- 



iment (a rat model for humans) to understand hypothermia therapy. Below we summarize 
the experiment as described by |Reynolds & Chiu| ( |2010| ). Thirty-eight approximately 8- 



week-old male Long-Evans rats were used in the experiment. Intraperitoneal transponders 
were implanted in the rats to record T c . Followed by a recovery time of 30 - 60 minutes, 
the rats were hemorrhaged from the carotid catheter with a constraint of a mean arterial 
pressure threshold of 40 mm Hg. The experiment was continued until the target shed 
blood volume (60% of the total blood volume) was achieved, which was then followed by 
a resuscitation intervention. Core temperature, T c , was logged by automated remote data 
collection every 15 seconds for the duration of the trial; T c for each rat was between 127 
and 246 time-steps (32 to 62 minutes) long. 



Reynolds & Chiu (2010) used broken sticks and bent cables to model the rat profiles, 
treating each as an individual time series. As such, they needed to omit some "outlying" 
profiles that did not obviously follow the shape of either the broken stick or the bent cable, 
and to truncate some other profiles which violated linearity of the incoming or outgoing 
phase; analyses for only 23 rat profiles were reported. In contrast, with our more general 
mixture methodology (Section[3]) proposed here, we can unify the inference from all 38 rats 
to address questions of broad interest about the underlying rat population, particularly: (a) 
How long did it take for the T c trend to show an obvious change because of hypothermia? 
(b) What were the rates of increase/decrease before and after the change? (c) What was 
the time point at which the trend went from increasing to decreasing, or vice versa? Like 



Reynolds & Chiu (2010), we consider data from the start of hemorrhage until resuscita- 
tion intervention. On the other hand, an obvious advantage of our method is that properly 
accounting for the longitudinal context allows pooling of information from the entire sam- 



pie, overcoming computational difficulties due to the apparent violation of the shape of the 
broken stick or bent cable for certain individuals. 

3. The Flexible Mixture Longitudinal Bent-Cable Model 

When m individuals can be regarded as having been randomly selected from some pop- 
ulation and repeated measurements are observed for each individual, it is useful to unify 
information from all m to aid the understanding of the population as well as subject-specific 
behaviour. For the i th individual (i = 1,2, ... , m), let there be rii measurements, and let 
tij denote the j th measurement occasion, j = 1,2,..., We model the corresponding 
response at time Uj, denoted by y^, by the relationship 

VH = f(fij, 0%) + ey (1) 

where is the vector of regression coefficients for the ith individual, /(■) is a function 
of Uj and 6i to characterize the trend of the subject-specific data, and e# represents the 
random error component, which accounts for measurement error and possibly additional 
within-individual error. 

For the types of data under consideration, in light of the apparent three phases - linear 
incoming and outgoing, and the adjoining curved transition - we characterize the individual 



profiles by the bent-cable function (Chiu et al. 2006), given by 



f(tij, 0i) = (3 i + PuUj + f3 2i q(ti j ,cxi), (2) 



where °^ = a~. ~ Ti ' - 7 ^ + ^ ~ r i) 1 {% ~ T i> T<} ( 3 ) 



7~i Ti 

with (3 i = (fioi, pu, f3 2 i)' and a, = (74, Ti)' being the vectors of linear and transition co- 
efficients, respectively, and 6i = (/3-, a-)'. For each individual i, /3 0i and f3u are, respec- 
tively, the intercept and slope of the incoming phase; (3u + (3 2 i, the slope of the outgo- 
ing phase; and t% and ji, the transition parameters which represent the center and half- 
width of the bend, respectively. Henceforth, we will denote f(tij,&i) and q(tij,cti) by 
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fij and qij. Note that % = reduces the bent cable to a broken-stick model for which 
Qij = (*v _ r *) 1 {^i - > 0} (Figure gb)). 



The critical time point (CTP), as defined by Chiu & Lockhart (2010), is the time at 
which the slope of the bent cable changes signs (Figure [2]). Thus, for a gradual transition, 
the CTP is Tj — 7« — 2/9i i 7i//5 2 i- Note that this formula is not meaningful when the slope of 
the cable does not change signs. When ^ = 0, any sign change of the slope occurs at the 
point t ?; , the CTP for an abrupt transition. 

We consider a hierarchical mixed-effects modelling framework and regard 0i as ran- 
dom, through which we can obtain useful information regarding the questions: (1) How 
does the response change over time (a) individually and (b) at the population level? and (2) 
Do different individuals experience different patterns of change? Question (l)(a) charac- 
terizes each individual's pattern of change over time (commonly called within-individual 
or Level 1 variation), and (2) addresses the association between patterns of change (com- 
monly called between-individual or Level 2 variation). Additionally, there is a third level 
for Bayesian inference, which quantifies prior knowledge for (1) and (2). 

The framework as described thus far constitutes the longitudinal bent-cable model of 



Khan et aL\ ( |2009[ ). For a mixed-effects model, it is parsimonious in the sense that the 
underlying population model is the bent cable which has only five fixed-effects regression 
coefficients. However, undesirable shrinkage issues as described in Section [T] are evident 



when their framework is directly applied to the rat data (see Section 4.1 ). Thus, to avoid 
estimation bias due to shrinkage, we further assume that 

Al. each individual i potentially comes from one of two populations: Population A for 
which 7; = and Population G for which ^ > 0; and 

A2. each individual has probability to to have come from Population G (and, hence, prob- 
ability 1 — ijj from Population A). 
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Conditional on the random effects 0j's, the sets of repeated measurements {yn,yi2, • • • , 
y irii } and {yki, yk2, ■ ■ ■ , ykn k } are assumed independent for % ^ k. To account for additional 
serial correlation among y^'s remaining after what has been accounted for by the Oi's, we 
assume at Level 1 that e^'s follow a stationary AR(p) process with a common p: 

(4) 



h£i,j-l + <P2ti,j-2 + ...+<Pp€ 



•i ,y—p Vij i 



where (p = 2 , ■ • • > 0p)' is the vector of AR(p) parameters, and [%-|of] ~ iV(0, of) 
for all j = 1, 2, ... ,n j. Furthermore, we consider a conditional likelihood framework for 
Level 1, where the initial p observations for each i, yj 1 ^ = (yn,yi2, yi P )', are treated as 
known, whereas y- 2 - 1 = (y ijP+ i,y iiP+ 2, ■ ■ ■ , y^m)' are random. This framework for y^ 1 ' and 
y f was also considered by 



Chiu & Lockhart 



(2010) for frequentist bent-cable regression 



for a single profile, and by Chib ( 1993 1 in a Bayesian approach for linear regression. 

Assumptions Al and A2, together with Equations ([l])-((4]), constitute our flexible mixture 

longitudinal bent-cable model. Letting: 

v p 

t%j ^ 4*k t%,j—ki Tij Qij ^ 4*k Qi,j—ki 



k=l 
P 



k=l 



and 



IMi = /3pj(l -} J <t>k) + PuXij + faiTij + <pk Vi,j-k 



k=l k=l 

for j = p + 1 , p + 2, . . . , rii, our choices of distributions for the relevant quantities allow us 
to rewrite the model as 

[yi,p+t\yi,t, Vi,t+i, Vi, P -i+t, ®h <t>, of] ~ N(ni jP+u ai)Vt=l,...,rii-p, (5) 

g(a l \I i ) = (l-I i )LN(r l \ii TA ,al A ) + I l LN 2 (oi i \^ a ^ a ), I, (6) 

Ii ~ BER(u) J 
[M^|hi,H!] ~ iV 8 (hi,Hi), [M a |h 2 ,H 2 ] ~ iV 2 (h 2 ,H 2 ), 
[0|h 3 ,e 3 ] ~ Ap(h 3 ,e 3 ), | ocflx] ~iV(a ,ai), 
[E^>x,Ax] -^(//^(z/xAi)- 1 ), [S-> 2 ,A 2 ]~ W^,^)- 1 ), J>. (7) 

[a" 2 |6 ,&i] [a- 2 |do,rfi] ~ G( l - 

[a;|c ,cx] ~ S(c ,cx) 



'do dj\ 
2 ' 2 /' 
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where fi g = (/i , ^1,^2)' and are, respectively, the mean and covariance of (3^ fi TA and 
a 2 TA are the mean and variance of log (tj) for Population A, fi a = (/z 7 , /x T )' and S a are the 
mean and covariance of log (cti) for Population G; and N p , LN P , BER, W, G and B stand 
for p-variate normal, p-variate lognormal, Bernoulli, Wishart, gamma and beta distribu- 
tions, respectively. Levels 1 and 2 are ([5]) and ([6]), and Level 3 is |7]) with the hyperparam- 
eters assumed known (see supp. material Section SI). Note that the distribution of ctj is a 
mixture of a univariate and a bivariate lognormal distribution corresponding to assumption 
Al; a.i = (74, Ti)' for Ii = 1, and cti = [tj] for U = due to a deterministic 7^ = 0. 

4. Rat Data Analysis: Bayesian Inference and 

Implementation 

4.1 The Longitudinal Bent-Cable Model 

The assumption that the samples come from Population A only is, perceivably, a restrictive 
and unrealistic assumption for a physiological phenomenon. As the existing framework 



by Khan et al. (2009) (Model G) allows an arbitrarily small 73 > for each i, we first 
applied it to our rat data to generalize this restrictive assumption. We observed an unusually 
large upper limit for the 95% credible interval for (E a ) n , i.e., the variance of 7$. This 
impracticality can be explained by noting that the presence of any rat i whose posterior 
draws for 7$ are arbitrarily small (e.g., < 1CT 3 ) can substantially inflate the corresponding 
posterior draws for (T, a ) n . As several rats exhibit a virtually abrupt transition while others 
do not, this resulted in an unreasonable estimate of (E a )n. 

4.2 The Flexible Mixture Longitudinal Bent-Cable Model 

Thus, there is practical need to generalize Model G by further extending it to a mixture of A 
and G as described in Section[3j As we explain below, the analysis using our mixture model 
provides strong evidence that supports the existence of not just A nor just G for the rat study. 
Bayesian inference is carried out by Markov chain Monte Carlo (MCMC), where we 
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sample from the posterior distribution by the Metropolis within Gibbs algorithm (Smith & 



Roberts, 1993). For implementation, we work out the full conditional for each parameter 
(see supp. material Section S2). We employ the Metropolis algorithm to draw samples of 
CKj, the sole parameter for which the full conditional can be expressed only up to a propor- 
tionality constant. The full conditional for (p is Gaussian; we take the proportion of draws 
(from the full conditional for (p) for which stationarity is satisfied as an estimate of the con- 



ditional probability of stationarity for the AR process (Chib 1993). We consider several 
models assuming {e^ } to be AR(p) for p = 0, 1, . . ., and choose the one for which the esti- 
mate of the deviance information criterion (DIC) is minimum ( |Spiegelhalter et al. 2002 ). 



Since our assumption for the ot^s involves lognormal distributions, we can use Level 2 
medians, namely M 7 = exp {^ 7 } and M T = exp {fi T } for Population G and M TA = exp {/jl Ta } 
for Population A, to describe the transition locations. We can also use Level 2 standard 



deviations of ji and 7* for G, namely <S 7 = y/exp{2pj + (E a )ii} x [exp {(£ a )n} - 1] and 
S T = \/exp{2/i T + (S a )22} x [exp{(S a ) 2 2} - 1] to describe the between-individual variabil- 
ity of these transition parameters. Posterior means or medians of Ais and Ss can be easily 
approximated using the MCMC samples. 

We proceed to analyze the aforementioned rat data using our flexible mixture bent- 
cable approach. We denote time by t^, j = 1,2, .. . , n^, where tn = refers to the 
starting point of the study for rat i (i = 1, 2, . . . , 38), and each subsequent time increment 
is 15 seconds. Any parameter estimate (Level 1 or 2) is based on the posterior mean or 
median, depending on the extent of asymmetry of the corresponding marginal posterior 
density. Note that 61 has its own posterior distribution, inducing a posterior distribution 
for the bent cable at each observed t^. Therefore, we regard the MCMC sample mean 
of fij as the fitted value fy. Individual-specific fitted curves are then interpolated based 



on the values; see Appendix A.l A fitted population curve is produced based on the 



estimates of the theoretical medians for ^ and a, : from Level 2. Similarly, we define the 
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CTP for Population G as M T — M 7 — 2^\M 1 / H2', thus, we use the posterior mean of this 
expression to make inference for this CTP. Estimates for the other parameters for Level 2 
Population A/G theoretical medians/standard deviations are produced similarly. 

5. Results 

Initially, we consider several flexible bent-cable models based on the degree of within- 
individual dependency among the repeated measurements, which is assumed through an 
AR(p) process for p = 0, 1, 2, 3. Model selection procedure reveals a smallest DIC for the 
AR(0) assumption. Fixing p = 0, we then analyzed the data using Models G and A (i.e., as- 
suming that the sample arises from Population G only and Population A only), and observed 



the smallest DIC for the proposed flexible model; see Appendix A. 2 for details. Therefore, 
we report here the results for the flexible model with AR(0) within-individual noise. 

Some posterior characteristics of parameters for the two populations are given in Ta- 
ble [T] and the population fitted curves are displayed in Figure [3] The posterior mean 
for u) is 0.39, suggesting that about 39% of the rats belong to Population G which exhibits a 
gradual change in T c . Posterior means for M T ±M y are 10.11 and 29.03 minutes, implying 
that the population transition begins approximately 10.11 minutes from the time of hem- 
orrhage and lasts for about 18.92 minutes, followed by a significant linear decrease at the 
rate of 0.013°C per 15 seconds (the posterior mean for /ii + /x 2 is —0.013 with 95% cred- 
ible interval (—0.016, —0.008)). The remaining 61% of the rats, approximately, exhibit an 
abrupt linear decrease at the same rate from the transition time point. We also see a signifi- 
cant linear increase in population T c at the rate of 0.003°C per 15 seconds in the incoming 
phase (95% credible interval of the incoming slope is (0.001, 0.006)). Moreover, virtually 
identical metabolic thresholds associated with a breakdown in the compensatory mecha- 
nisms for the two populations are observed (see Figure [3]): posterior means for Population 
G and A CTPs are 14.28 and 13.89 minutes, respectively (Table [T]). Thus, for G, the drop 
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in T c started at approximately 14.28 minutes after hemorrhage, and 13.89 minutes for A. 

[Table 1 about here.] 
[Figure 3 about here.] 

Figure [T] shows examples of the individual fitted curves. For all but one of the six rats 
displayed (and all but a total of four out of the entire sample of 38 rats), the fits appear 
very reasonable as the observed data closely agree with the respective fitted lines, and 
the estimated transitions ff and r ± 7) demonstrate that our methodology picks up the 
two types of transition adequately. The remaining rat (Figure [TJc)) appears to be unusual, 
exhibiting linearly decreasing trends throughout (recall Section[T]); again, this is one of four 
rats among the 38 who do not cleanly fall into either population. With our methodology, 
these four are estimated to have arisen from Population A. Given our small dataset, we do 
not consider a potential third population to avoid overfitting. 

The posterior characteristics of the theoretical standard deviations and correlations are 
given in Table [2j Since the biological conditions of different rats should vary to some ex- 
tent, we can expect some variation in T c 's at the time of administering hemorrhage. This 
is reflected in the estimate of the standard deviation of which is 0.535. After adminis- 
tering hemorrhage, we see very little variation in the slope parameters (estimated standard 
deviations for f3u and fin + fai are 0.008 and 0.011, respectively), that is, all rats exhibit 
very similar rates of increase/decrease before/after the transition period. Significant nega- 
tive correlation between (3 U and (3 2 i (corr((3u, f3 2 i) = —0.476 with 95% credible interval 
(—0.711, —0.204) which excludes zero) indicates that the steeper the incoming slope going 
up, the bigger the drop in slope for the outgoing phase. 

[Table 2 about here.] 
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From Table |2} we see considerable variability in the times to maximal T c and the vari- 
ability in the times to transition zones. This fact is reflected in the posterior medians for the 
standard deviations of 7$ and 7* for Population G (7.643 and 10.172 minutes, respectively), 
and of n for Population A (9.312 minutes). We also see corr^i, 7$ — 7$) = —0.816 with 
95% credible interval (—0.976, —0.602), which excludes zero. Significant negative corre- 
lation between 7^ and — 7$ indicates that for individuals from Population G, the sooner the 
gradual transition takes place, the wider the transition zone so that there will be a delayed 
linear drop in the outgoing phase, and vice versa. 

In summary, our analysis yields the following points of clinical interest: (a) about 61% of 
the rats exhibit an abrupt linear drop in T c during hemorrhage, whereas the remaining 39%, 
approximately, exhibit a gradual transition followed by a linear drop; (b) all rats are from 
populations that show approximately the same metabolic threshold (about 14 minutes after 
hemorrhage) associated with a breakdown in the compensatory mechanisms; (c) either pop- 
ulation shows a significant increase of T c followed by a significant decrease; (d) all the rats 
exhibit very similar rates of increase and decrease in T c before and after the transition pe- 
riod, respectively; (e) there is a considerable amount of between-rat variability in the times 
to maximal T c and transition zones; (f) the sooner the gradual transition takes place, the 
wider the transition zone, and vice versa; (g) although assuming within- subject conditional 
independence may be unrealistic for some problems, we demonstrate in Scenarios 2 and 3 
in the next section that points (a)-(f) above should be reasonably robust to this assumption. 

6. Simulations 



First, we supplement the motivation for our mixture methodology as seen in Section 4.1 



That is, we show the importance of hypothesizing Population A in addition to G in a more 
general context, despite that the abrupt broken stick is the limiting case of the smooth bent 
cable. To this end, we present Scenario 1, where we fit Model G when, in reality, both 
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Populations A and G exist, with G heavily dominating A: (a) tu = 0.90 and (b) u = 0.95. 
In both la and lb, {e^} is an AR(1) process with = 0.70, where p = 1 is treated as 
known when fitting Model G. 

Second, to illustrate the fact that our flexible methodology can perform well with re- 
spect to the population regression coefficients even for a misspecified correlation structure 
for {ejj}, we present Scenarios 2 and 3, where {e^} has AR(1) or AR(2) structure. In each 
case, we analyze the data assuming p = 0,1, and 2, and that the samples come from two 
potential populations (A and G). In all the scenarios, we take m = 20, n = rii = 150 for 
% = 1,2, . . . ,m and tij — j — 1 for j = 1,2, ... ,n. Model parameter values were chosen to 
allow reasonable generalization, and are given in Tables [3]- [8] in supp. material Section S|4| 

For each simulation, 500 data sets are generated, and 100,000 MCMC iterations are 
used to approximate posterior distributions per set. Posterior summaries are averaged over 
the 500 sets for each parameter, and the coverage probability of 95% credible intervals 
(proportion of such credible intervals out of 500 that capture the truth) is calculated. 

6.1 Results for Scenario 1 

Numerical results are tabulated in Tables [3] -[5] in supp. material Section S4.1. We see that 
Model G performs well with respect to all but one population regression coefficient, /i 7 , 
the bend's half-width (for Population G). Specifically, the average of the posterior means 
for each parameter except /i 7 is close to the true parameter value, and the corresponding 
coverage probabilities are all reasonably close to the nominal 0.95. When u = 0.90, we 
see underestimation and under coverage for /i 7 . This can be explained by noting that the 
average of the posterior means for each 7, is expected to be approximately zero for profiles 
that originate from A; this leads to underestimation of the population counterpart Note 
that if we would model this data set using our flexible methodology, /i 7 would represent 
only the profiles that originate from G, and hence, underestimation for /i 7 would not be 
expected, and coverage for /x 7 would be close to the nominal 0.95. Indeed, this is evident 
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from the results for tu = 0.95 (fewer abrupt profiles than tu = 0.90): the average of the 
posterior means for /x 7 is 50% closer to the true value of 3, and also the coverage is 67% 
closer to the nominal 0.95 (see supp. material Section S4.1). 

Details about E^, E a , and of 's also appear in Section S4.1. The main conclusion is that 

(i) the estimation of (£ a )n and (E a ) 2 2 is more accurate for u « 1, and (ii) misspecifying 
the model as Model G, when in reality both populations A and G exist, negligibly affect 
the estimates of erf's or E^. 

The above simulation results demonstrate the importance of modelling Population A 
distinctly from G using our flexible methodology to analyze data that resemble those from 
the rat experiment; note the rat experiment required an even more extreme need for a mix- 
ture, with 95% credible interval for cu being (0.23, 0.55). 

6.2 Results for Scenarios 2 and 3 

Numerical results for Scenarios 2 and 3 are given in Tables [6]- [8] in supp. material Section 
S4.2. Our methodology performs well for both scenarios with respect to the population 
characteristics: averages of posterior means are all close to the true parameter values, and 
coverage probabilities (from 0.92 to 0.99) are all reasonably close to the nominal 0.95. This 
suggests that our Bayesian inference for population characteristics is robust to ignoring 
certain types of serial correlation. 

Details about E^, E a , <r\ A , and of 's also appear Section S4.2. The main conclusion is 
that (i) the estimation of E^, E Q , and a\ is quite accurate for correctly specified models, 
though underspecifying p may result in under coverage for (E Q )n, (E Q ) 22 and a* , and 

(ii) an underspecified p leads to overestimation of of. In particular, we observe very poor 
coverage for of if we incorrectly analyze a data set by an AR(0) assumption when, in 
reality, it exhibits serial correlation over time. Although poor coverage may not be ideal 
in certain cases, of primary practical concern in the rat analysis is the inference for the 



17 



population regression coefficients, for which our methodology demonstrates robustness. 

7. Conclusion 

Induced hypothermia potentially saves lives under physiological trauma. Yet, without ex- 
treme care, it can also threaten survival. Therefore, controlled administration of hypother- 
mia is of paramount importance. In this article, we developed the flexible mixture bent- 
cable framework to quantify the transition of core body temperature, T c , during induced 
hypothermia in a rat model for humans. Our analysis reveals important clinical informa- 
tion that can be very valuable in administering hypothermia therapy. Aside from crucial 
information at the population level, another aspect of our longitudinal framework which 
clinicians would find valuable is the inference for the temporal trend exhibited by individ- 
ual observational units: the current inference for an individual may provide guidelines to 
future administration of hypothermia therapy to the same individual. 

The most appealing feature of our method may be its greatly interpretable parameters, 
and that useful information can be obtained at the small cost of estimating very few fixed- 
effects regression coefficients. Moreover, pooling information from many individuals leads 
to shrinkage, so that mild deviations of observed profiles from the broken-stick/bent-cable 
structure do not hinder model fitting; in contrast, deviations considered mild can render 



the single-profile bent-cable regression method infeasible (e.g., |Reynolds & Chiu[ |2010). 
Despite the broken stick being the limiting case of the bent cable, reliable inference for 
the underlying population transition of T c requires that the stick population be an explicit 
component of a mixture model comprising both stick and cable populations, even if the 
cable population dominates in size. Moreover, the mixture allows better inference for the 
CTPs for separate populations (not presumed identical a priori); we have evidence that the 
population in the rat study consists of more than just A or just G, so that the inference (for 
the CTP and other parameters) would be incorrect if we did not use the mixture. There- 
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fore, our extension of bent-cable regression to model longitudinal data for multiple units 
provides a desirable statistical tool to characterize a special type of continuous temporal 
trend — one showing a change due to a shock that exhibits both gradual and abrupt tran- 
sitions. Although it would require further subject-matter research to investigate the phys- 
iological reason for certain individuals to exhibit an abrupt instead of gradual change, our 
flexible bent-cable approach offers an empirical solution for identifying them and making 
integrated inference for their CTP alongside individuals who exhibit a gradual change. Our 
methodology, under a general regression modelling framework, can classify observational 
units in the same longitudinal study as exhibiting either an abrupt or gradual transition. 
It provides not only inference that is more realistic, but also insights into the underlying 
behaviour within a population. As such, it is applicable to the rat model for induced hy- 
pothermia, and potentially to a wide variety of other situations. Also, if there were enough 
observations to support, say, a third population, our method could be easily extended to 
include a third component of the mixture model. 

Our method is intended for only stationary AR(jo), p > 0, processes for {ey}, though 
simulations suggested that assuming an AR(0) structure even when serial correlation ex- 
ists among repeated measurements does not lead to problematic bias when characterizing 
the populations. In this case, serial correlation is induced by the random regression coeffi- 
cients. Some directions of extension to address this and other limitations are suggested in 



Appendix A. 3 Overall, the flexible mixture bent-cable model for longitudinal data as pro- 
posed in this paper has many attractive properties and has allowed us to model data from, 
and provide informative interpretations for, a scientific problem of great practical interest. 
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Appendices 
A.l Approximating Fitted Values 

The parameter vectors /3 i and aij have their own posterior distributions, so the bent-cable 
function fa itself has a posterior distribution at each observed time point t^, j = 1, 2, . . . , n,. 
We consider the posterior of the bent-cable function to produce the fitted values by taking 
the MCMC sample means of the bent-cable function. So, the bent cable for the i th individ- 
ual at observed time is = (3 0i + flu + 02%Qij^ an d the corresponding fitted values are 

T 

T 



with gg 



s=l 

where T is the length of the MCMC samples. 

A.2 Model Selection 

Model selection procedure is carried out by comparing DICs. We initially consider four 
flexible mixture bent-cable models assuming e^-'s to follow an AR(p) process for p = 
0, 1, 2, 3. Note that we consider a conditional likelihood framework for an AR(p) process, 
where the initial p observations for each i are treated as known. Therefore, as suggested 
by |Chiu & Lockhart (2010), the analyses were initially performed on a reduced dataset to 



make the DICs comparable for p = 0,1, 2, 3. Specifically, we consider (y^, yi$, yi, ni )' 
as the response vector (random) for the i th individual for all comparisons, while (y^), 
(z/i,2) yi,s)' an d {yi,x,yi,2, 1/1,3)' are treated as known for p = 1, 2 and 3, respectively. That is, 
we dropped the first 3 — p observations for each i for p = 0, 1, 2, 3, respectively. 

Preliminary analysis (not shown) reveals that the data exhibit nonstationarity when as- 
suming p > for {eij}: the proportion of draws from the full conditional of for which the 
stationarity condition is satisfied is close to zero. For example, an AR(1) assumption with 
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prior ~ N(0, 10 4 ) leads to = 0.99 with DIC ^ 5.16 x 10 s . To achieve stationarity, we 
also consider ~ N(0, 5 x 10~ 5 ) and <p ~ JV(0, 2.5 x 1(T 5 ) that lead to DIC « 10722 
and —14004, respectively. In addition, the fitted coefficients change depending on the prior 
variance for 0. Nonstationarity was also observed for AR(2) and AR(3) assumptions. 

In light of the extreme sensitivity to the prior for <p while assuming stationarity, we 
assume within-individual conditional independence (AR(0)), such that within-individual 
dependence among repeated measurements is due solely to the inclusion of the random 
effects Oi's. Our simulations (Section[6]) reveal that though the estimates of the a/s could 
be less reliable, the flexible methodology can perform well with respect to the population 
parameters even for a misspecified correlation structure for the e^-'s. Since our main goal 
is to make inference about the populations, we report results for AR(0) with DlCpiexibie ~ 
—14729 (smallest observed) in Section|4j We also analyzed the data using Models G and A 
for AR(0), for which DIC G « -14570 and DIC A ~ -13819, that is, our flexible mixture 
bent-cable model yielded better goodness of fit. Finally, note that the reported inference is 
actually based on the full data, i.e., using the reduced dataset as described above was solely 
for the purpose of DIC comparisons. 

A.3 Possible Extensions 

Although tailored for the rat study, the mixture longitudinal bent cable is perceivably ap- 
plicable to other studies involving profiles that exhibit abrupt and/or gradual transitions of 
temporal trend. Thus, extensions of our framework may be desirable in some cases. For 
example, our framework is intended for only stationary AR(p), p > 0, processes for within- 
individual noise, and it might be useful to extend the framework to specifically account for 
nonstationarity. Other possible extensions include (i) with sufficient data, allowing for 
additional populations to be considered in the mixture, and (ii) allowing the variation of 
profiles to depend on both random and systematic components (covariates). 
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SUPPLEMENTARY MATERIAL 
SI. Choice of the Hyperparameters 

Values of the hyperparameters reflect our prior knowledge. When little is reliably known 
about the individual trajectories beyond its functional form of the bent-cable, it is rea- 
sonable to choose the hyperprior values that lead to fairly vague, minimally informative 
priors ( |Carlin[ |1995| ). 

The choice of a mean vector (e.g., hi, h 2 or h 3 ) has very little effect on Bayesian esti- 
mation, as long as the respective variance parameters (diagonal elements of Hi, H 2 or H 3 , 
respectively) are taken to be very large which lead to flat priors (Song] [2007 ). Therefore, a 
common practice is to choose a zero mean vector and a covariance matrix, say, Hi such that 
Hf 1 « O, where O is a matrix with all its elements zero (Davidian and Giltinan, 1995 1. 



We use the parameterization of the gamma distribution as given in Chib ( 1993 1. For 



example, [<7^ 2 |c?o>^i] ~ G( y , y). Small values of the hyperprior parameters (e.g., do = 
d\ = 10~ 4 ) lead to a diffuse prior. 

We use the parameterization of the Wishart distribution as given in Carlin (T995). For 
example, [El 1 ^, Aj] ~ W{y\, {y\K^)~ x ). Setting the degrees of freedom equal to the 
order of the scale matrix (e.g. 3 for the prior of Ea 1 ) makes a Wishart prior nearly flat 



(Wakefield et ah, 1994). The matrix Ai (or A 2 ) is chosen to be an approximate prior 
estimate of E^ (or E a ). In the absence of such prior knowledge, one may use the sample 
covariance matrix of the individual- specific estimates of the regression coefficients; the R 
( |R Development Core Team} [2011) library "bentcableAR" (jChiu} |2008-2010[ ) for single 
profile bent-cable regression can be useful in this regard. 

Since < to < 1, we choose the beta distribution [<jj\cq, ci] ~ B(cq, ci) in our model. 
In the absence of prior information, one may choose c = c\ = 1 which leads to U(0, 1) 
distribution. 
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S2. Full Conditionals 



For the full conditionals of the model parameters, let 




Qij ~ Efc=l 0k Qi,j-k, 

(r ijP+1 ,.. .,r i>n .)' and 



Xi = (1 - ELi^fe, Xj, ri); 

• e u = Uij ~ Poi ~ Pa t ij ~ Pn Qij for 3 = V + 1, • • • , ^ = • • • , Q,nJ' and 
V -1 = E™ i °i~ 2 w 'i w * + H 3 where w * is a{ni-p)xp matrix with the k th row 
given by (e ijk +p-i, u,k+p-2, • • • , Q,fc); 

• m A = YZi (! - J i) and m c = E™ i 

• ^ = logai = (log7i,logri)' and = logr^; 

• 3 = E™i ft. I = Er=i * 6, and « = Er=i (i - 4) 

. M- 1 = of X$ X, + E^ 1 ; 

• U^ 1 = m S, 1 + H^ 1 and U 2 1 = m G S" 1 + H 2 \ 



An appealing feature of the bent-cable function is that it is partially linear - given a i5 
f(Uj, 0i) is linear - and we can exploit this fact to derive a closed-form full conditional 
for /3j. However, the full conditional of oli can be expressed only up to a proportionality 
constant, and is given by 



7r(a;|.) oc exp I - ^(z* - X; ft/fo - Xj ft) j x ^- exp | - - aO 2 





The full conditionals of the remaining parameters are 
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o T l n + a 1 l a 1 
v m A cr rA 2 + a x m A a T ^ + a x ! 



[fi a \.\ ~ iv 2 (u 2 (s- 1 1 + H 2 1 h 2 ), U 2 ^ 
m G + ^ 2 , [X)/i(^-A* a )(^-A* a )' + ^A 2 ] J 



oZ 2 \.]^G 



i=i 



-., i-j \ 2 ' 2 

rii-p + do (zi - Xf /3j (zi - X; f ) + d x 



i=i 



~ 5(m G + c , m A + ci). 

S3. Fitted Population Curves for the Rat Profiles 

In the main text, we presented the rat data analysis; the fitted population curves were dis- 
played in Figure 3. Here, in Figure [4j we reproduce the same figure but overlaid with all 38 
rat profiles; the population fitted curves are displayed in bold. Figure [4] displays the whole 
range of shapes of the rat profiles. It also shows that the profiles are well represented by 
the population fitted curves. 

[Figure 4 about here.] 

S4. Detailed Simulation Results 

Main findings of our simulations were summarized in the main text. Here we present the 
numerical results with supplementary information. 
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S4.1 Scenario 1 

Simulation results for Scenario 1 are presented in Tables [3] - [5]. In the main text, we dis- 
cussed our findings for the population regression coefficients of Table |3j 

[Table 3 about here.] 
[Table 4 about here.] 
[Table 5 about here.] 

We calculate the percentage closer to the true value for /i 7 (given on page 16 of the 
main text) as follows. The true /i 7 is 3.00, whereas the averages of the posterior means are 
2.88 and 2.94 for to = 0.90 and 0.95, respectively. Then, the average of the posterior means 
when u = 0.95 is 100{(3.00 - 2.88) - (3.00 - 2.94)}/(3.00 - 2.88) = 50% closer to the 
true value when compared to = 0.90. Similarly, the coverage for /i 7 is 100{(0.95 — 0.86) — 
(0.95 - 0.92)}/(0.95 - 0.86) « 67% closer to the nominal 0.95. 

From Table [4j we see that coverage probabilities for the elements of are all close to 
0.99. Since T* a takes into account both A and G, large variabilities among the 7/s and r/s 
are expected. In our simulation study, large variabilities are indeed reflected through the 
overestimation for each of the variance parameters (£ a )n and (T> a ) 2 2- Because of the over- 
estimation, we also see low coverage probabilities. Moreover, for a true (£ a )ii of 0.020, 
the average of its posterior medians is 0.114 when cu = 0.90, and is 0.067 when to = 0.95. 
Hence, estimation of (£ a )n is more accurate for to close to 1; the same is true for (S a ) 2 2- 

Finally, the average of the posterior medians for each o\ is close to the truth, and cover- 
age probabilities are all close the the nominal 0.95; see Table [5j That is, misspecifying the 
model as Model G, when in reality both populations A and G exist with G being dominant, 
has virtually no effects on the estimates of erf's. 
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S4.2 Scenarios 2 and 3 

Simulation results for Scenarios 2 and 3 are presented in Tables [6] - [8j We discussed our 
findings for the population regression coefficients (Table [6]) in the main text. 

[Table 6 about here.] 
[Table 7 about here.] 
[Table 8 about here.] 

For S Q and of , coverage probabilities are all close to 0.99 (slight over coverage) 
for correctly specified models (Table|7]). However, we observe poor coverage for (£ a )n for 
data sets generated from AR(1) and AR(2), but using an AR(0) fit: coverage probabilities 
are 0.78 and 0.80, respectively. For such model mis specification, we also see a slight 
under coverage for (T, a ) 2 2 and overestimation of (S Q )n, (S a )22 and a* This suggests 
that underspecifying p as zero may result in overestimation of transition parameter prior 
variances. Over-coverage, as we have observed for some parameters in all three scenarios, 
is of much less concern in practice than under-coverage. 

Finally, we see noticeable differences in the estimates (average of the posterior medi- 
ans) of of 's for different p's. In general, an underspecified p leads to overestimation of 
of. We observe very poor coverage (from 0.02 to 0.04) for of if we incorrectly analyze a 
data set by an AR(0) assumption when, in reality, it exhibits serial correlation over time. 
However, the problem is much less severe for an underspecified p that is positive. Al- 
though such poor coverage may not be ideal in certain cases, of primary practical concern 
in the rat analysis is the inference for the population regression coefficients, for which our 
methodology demonstrates robustness. 
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Table 1: Posterior summaries for the two populations of rats assuming AR(0) noise: posterior means 
for the population slope parameters (fii and ^2) are in "per 15 seconds" and those for the population 
transitions are in minutes. 





Posterior 
mean 


95% credible 
interval 


uj (Probability of being from G) 


n on 
I). 39 


(U.23, U.55) 


Ho (Incoming intercept) 


67. So 


(37 .21, 37.5b) 


[i\ (Incoming slope) 


U.UU3 


(U.UU1, U.UUb ) 


Li '2, \^ lAAt/lt/H^t/ ULLWLL11 clllLl \JLl lii^Jlllii SHJJJ^Sy 


— n m fi 


C_n n2fl -d (1111 


M TA (Population CTP for A) 


13.89 


(10.59,17.34) 


7W 7 (Half-width of the bend for G) 


9.46 


(5.45,13.62) 


M T (Center of the bend for G) 


19.57 


(13.67,25.30) 


M. r ± M. 1 (Transition period for G) 


10.11 to 29.03 




Mr - M 1 - 2^iA^ 7 / / u 2 (Population CTP for G) 


14.28 


(6.33,21.84) 
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Table 2: Rat data analysis - posterior summaries of the standard deviations and correlations asso- 
ciated with Sy3, T> a and a^. ; posterior medians for the standard deviations of the linear parameters 
(Pu and /?2i) are in "per 15 seconds" and those for the transition parameters (7$ and are in 
minutes. 





Posterior 
median 


95% credible 
interval 


S.D. of fa 


0.535 


(0.423,0.669) 


S.D. o£0u 


0.008 


(0.006,0.010) 


S.D. of/3 2i 


0.123 


(0.010,0.016) 


S.D. olpu + fa 


0.011 


(0.009,0.014) 


Corr. between fa and /3h 


0.023 


(-0.296,0.343) 


Corr. between fa and /?2i 


-0.001 


(-0.323,0.319) 


Corr. between j3u and /?2i 


-0.476 


(-0.711,-0.204) 


S.D. of 7j for Population G 


7.643 


(2.967, 19.068) 


S.D. of Tj for Population G 


10.172 


(4.728,20.928) 


Corr. between 73 and Tj — 7$ for Population G 


-0.815 


(-0.976,-0.602) 


S.D. of Ti for Population A 


9.312 


(5.030,16.200) 
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Table 3: Simulation scenario 1 results with n, = 150 for all i and m — 20: average of 500 posterior means 
of the population regression coefficients and the AR parameters; also coverage of 95% credible intervals. 







Simulated u = 0.90 


Simulated ui — 0.95 






Analysis using Model G 


Analysis using Model G 




True 


Mean, Coverage 


Mean, Coverage 


Mo 


244.00 


244.33,0.95 


244.51,0.96 


Mi 


0.50 


0.49,0.93 


0.49,0.94 




-0.75 


-0.78,0.91 


-0.78,0.92 


H 


3.00 


2.88,0.86 


2.94,0.92 




4.00 


4.04,0.92 


4.02,0.95 




4.50 









0.70 


0.71,0.93 


0.71,0.93 
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Table 4: Simulation scenario 1 results with = 150 for all i and m — 20: average of 500 posterior means 
(medians for the variance parameters) of the variances and covariances (S^ and E a ) in the priors for the 
random regression coefficients; also coverage of 95% credible intervals. 







Simulated lj = 0.90 


Simulated = 0.95 






Analysis using Model G 


Analysis using Model G 




True 


Mean, Coverage 


Mean, Coverage 


(£/»)n 


125.00 


123.65,0.99 


123.38,0.98 


(5^/3)22 


0.03 


0.03,0.98 


0.03,0.97 


1 v ^ 

0/3J33 


U.Uo 


n HQ n no 

u.uo, u.yy 


n rtQ rt on 
U.Uo, U.»9 


( S ^)l2 


-1.00 


-0.95,0.96 


-0.94,0.97 


( S ^)l3 


0.50 


0.57,0.99 


0.57,0.99 


(5^/3)23 


-0.01 


-0.01,0.99 


-0.01,0.99 


(S Q )ii 


0.020 


0.114,0.69 


0.067,0.82 


(Sq,) 2 2 


0.030 


0.054,0.62 


0.043,0.78 


(^a)l2 


0.005 


-0.051,0.68 


-0.024,0.83 




0.050 
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Table 5: Simulation scenario 1 results with rij = 150 for all i and m = 20: average of 500 posterior medians 
of the innovation variances; also coverage of 95% credible intervals. 







Simulated us = 0.90 


Simulated cj — 0.95 






Analysis using Model G 


Analysis using Model G 




True 


Mean, Coverage 


Mean, Coverage 




0.34 


0.35,0.97 


0.35,0.94 




1.12 


1.14,0.95 


1.12,0.94 




1.75 


1.78,0.95 


1.76,0.96 


°l 


0.42 


0.42,0.95 


0.42,0.96 


°l 


0.74 


0.76,0.94 


0.74,0.94 


-I 


2.06 


2.08,0.95 


2.08,0.95 


a? 


1.16 


1.16,0.94 


1.16,0.93 


4 


1.28 


1.29,0.93 


1.27,0.93 


-I 


0.16 


0.16,0.95 


0.16,0.96 




0.77 


0.78,0.96 


0.77,0.94 




0.04 


0.04,0.96 


0.04,0.95 




0.03 


0.03,0.96 


0.03,0.94 


2 

<7l3 


0.91 


0.92,0.96 


0.92,0.95 




1.96 


1.97,0.94 


1.96,0.95 


*h 


0.32 


0.33,0.96 


0.32,0.96 




2.02 


2.02,0.95 


2.05,0.95 


°b 


0.89 


0.90,0.95 


0.90,0.96 




0.90 


0.90,0.94 


0.91,0.95 




0.82 


0.83,0.95 


0.83,0.95 


° 2 o 


2.89 


2.93,0.97 


2.91,0.96 
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Table 6: Simulation scenarios 2 and 3 results with = 150 for all i and m = 20: average of 500 posterior 
means of the mixing proportion ui, population regression coefficients and the AR parameters; also coverage 
of 95% credible intervals. 







Simulated e^ 's: AR(1) 


Simulated e^ 's: AR(2) 






Analysis assuming 


Analysis assuming 






AR(2) 


AR(1) 


AR(0) 


AR(2) 


AR(1) 


AR(0) 






Mean, 


Mean, 


Mean, 


Mean, 


Mean, 


Mean, 




True 


Coverage 


Coverage 


Coverage 


Coverage 


Coverage 


Coverage 


UJ 


A C A 

U.5U 


0.52, 0.98 


0.52,0.99 


0.52,0.97 


0.52, 0.97 


0.51, 0.95 


0.52,0.96 


Mo 


O A A AA 

244.00 


243.94,0.97 


244.45,0.94 


244.37,0.96 


244.64,0.96 


244.30,0.95 


244.45,0.96 


Mi 




0.50,0.97 


0.48,0.94 


0.49,0.95 


0.48,0.94 


0.48,0.95 


0.48,0.93 


M2 


-0.75 


-0.75,0.95 


-0.77,0.92 


-0.78,0.92 


-0.77,0.95 


-0.78,0.93 


-0.77,0.92 


M 7 


3.00 


2.95,0.95 


2.94,0.95 


2.95,0.93 


2.97,0.96 


2.95,0.97 


2.96,0.93 


Mr 


4.00 


4.02,0.98 


4.02,0.97 


4.04,0.96 


4.01,0.99 


4.02,0.96 


4.04,0.92 


Mr A 


4.50 


4.50,0.94 


4.50,0.96 


4.47,0.94 


4.49,0.95 


4.49,0.97 


4.47,0.94 


AR(1) <j> 


0.70 




0.71,0.93 






0.74,- 




AR(2) 0! 


0.80 


0.70,0.94 






0.80,0.95 






AR(2) </> 2 


-0.10 


0.005,0.92 






-0.10,0.95 
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Table 7: Simulation scenarios 2 and 3 results with n l = 150 for all i and m = 20: average of 500 posterior 
means (medians for the variance parameters) of the variances and covariances (cr%. , £^ and S Q ) in the priors 
for the random regression coefficients; also coverage of 95% credible intervals. 



Simulated e;/s: AR(1) Simulated e;/s: AR(2) 

Analysis assuming Analysis assuming 







AR(2) 


AR(1) 


AR(0) 


AR(2) 


AR(1) 


AR(0) 






Mean, 


Mean, 


Mean, 


Mean, 


Mean, 


Mean, 




True 


Coverage 


Coverage 


Coverage 


Coverage 


Coverage 


Coverage 




125.00 


126.13,0.98 


126.78,0.97 


124.32,0.98 


124.39,0.97 


123.11,0.98 


125.52,0.98 


(£9)22 


0.03 


0.03,0.99 


0.03,0.99 


0.03,0.98 


0.03,0.97 


0.03,0.98 


0.03,0.98 


(2,3)33 


0.03 


0.03,0.98 


0.03,0.99 


0.03,0.97 


0.03,0.98 


0.03,0.97 


0.03,0.99 


( S /3)l2 


-1.00 


-1.05,0.98 


-0.97,0.95 


-0.98,0.97 


-0.97,0.96 


-0.93,0.97 


-0.97,0.97 


( S /3)l3 


0.50 


0.49,0.99 


0.60,0.98 


0.53,0.99 


0.54,0.98 


0.53,0.99 


0.49,0.99 


(5^^)23 


-0.01 


-0.01,0.97 


-0.01,0.98 


-0.01,0.98 


-0.01,0.98 


-0.01,0.99 


-0.01,0.98 


(Sq,)h 


0.020 


0.021,1.00 


0.020,0.99 


0.059,0.78 


0.021,0.99 


0.036, 1.00 


0.058,0.80 


(Sa)22 


0.030 


0.031,0.99 


0.031,0.99 


0.045,0.91 


0.032,0.99 


0.032,0.99 


0.043,0.93 


(S a )i2 


0.005 


0.0002,1.00 


0.001,1.00 


-0.007,0.95 


0.001,0.99 


0.0003, 1.00 


-0.006,0.95 


< 


0.050 


0.57,0.97 


0.059,0.98 


0.069,0.98 


0.059,0.95 


0.059,0.97 


0.073,0.97 
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Table 8: Simulation scenarios 2 and 3 results with = 150 for all i and m = 20: average of 500 posterior 
medians of the innovation variances; also coverage of 95% credible intervals. 



Simulated e^-'s: AR(1) 
Analysis assuming 
AR(2) AR(1) AR(0) 



Simulated e^-'s: AR(2) 
Analysis assuming 
AR(2) AR(1) AR(0) 



True 



Mean, 
Coverage 



Mean, 
Coverage 



Mean, 
Coverage 



Mean, 
Coverage 



Mean, 
Coverage 



Mean, 
Coverage 



a\ 0.34 

a\ 1.12 

o\ 1.75 

a\ 0.42 

erf 0.74 

<t| 2.06 

a% 1.16 

o-| 1.28 

ct| 0.16 

a 2 10 0.77 

o\ x 0.04 
0.03 

a\ z 0.91 

of 4 1.96 

aj 5 0.32 

of 6 2.02 

cr^ 7 0.89 

of 8 0.90 

trfg 0.82 

4) 2.89 



cr 



12 



0.35,0.96 
1.12,0.94 
1.77,0.95 
0.41,0.94 
0.75,0.95 
2.08,0.96 
1.16,0.94 
1.29,0.95 
0.16,0.94 
0.78,0.97 
0.04,0.94 
0.03,0.95 
0.92,0.96 
1.99,0.94 
0.32,0.96 
2.03,0.97 
0.90,0.94 
0.90,0.94 
0.83,0.93 
2.92,0.93 



0.35,0.96 
1.13,0.97 
1.76,0.95 
0.42,0.95 
0.74,0.95 
2.08,0.94 
1.16,0.95 
1.30,0.96 
0.16,0.94 
0.79,0.95 
0.04,0.95 
0.03,0.96 
0.92,0.95 
1.95,0.95 
0.33,0.95 
2.03,0.94 
0.89,0.94 
0.90,0.94 
0.84,0.95 
2.92,0.94 



0.58 
1.89 
2.95 
0.71 
1.26 
3.52 
1.95 
2.17 
0.27 
1.32 
0.06 
0.06 
1.55 
3.36 
0.55 
3.40 
1.52 
1.53 



0.11 
0.10 
0.09 
0.08 
0.08 
0.09 
0.09 
0.08 
0.08 
0.08 
0.07 
0.09 
0.09 
0.06 
0.07 
0.07 
0.09 
0.09 



1.41,0.07 
4.86,0.10 



0.35,0.96 
1.12,0.94 
1.76,0.94 
0.42,0.96 
0.74,0.94 
2.09,0.95 
1.16,0.94 
1.28,0.95 
0.16,0.93 
0.78,0.95 
0.04,0.95 
0.03,0.94 
0.91,0.95 
1.96,0.96 
0.33,0.95 
2.04,0.96 
0.89,0.95 
0.91,0.96 
0.82,0.96 
2.93,0.96 



0.36,0.95 
1.16,0.95 
1.80,0.94 
0.43,0.92 
0.76,0.96 
2.10,0.94 
1.18,0.94 
1.30,0.96 
0.16,0.95 
0.79,0.95 
0.04,0.94 
0.03,0.95 
0.93,0.95 
1.99,0.94 
0.33,0.94 
2.05,0.95 
0.91,0.94 
0.92,0.96 
0.84,0.94 
2.99,0.95 



0.65,0.02 
2.09,0.02 
3.27,0.03 
0.78,0.02 
1.43,0.02 
3.86,0.03 
2.19,0.03 
2.43,0.03 
0.30,0.02 
1.44,0.03 
0.07,0.03 
0.07,0.02 
1.70,0.04 
3.62,0.04 
0.61,0.03 
3.85,0.03 
1.68,0.04 
1.70,0.04 
1.54,0.04 
5.47,0.03 
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(e) Rat 17 



37.0- 
36.5- 
36.0- 
35.5- 
35.0- 



38.4- 



38.2- 



37.8- 



(f) Rat 34 
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Figure 1: Observed profiles (grey curves) and the corresponding individual-specific fitted curves (solid, in 
black) along with 95% pointwise credible intervals (dotted curves) for selected rats; fits A and G virtually 
coincide for each of Rats 4, 9 and 17. Estimated transitions (i.e. f and r ± 7) are marked by solid vertical 
lines, and estimated CTPs (for Population G) by dotted vertical lines; the CTP estimate is not marked for 
Rat 3 because the estimated slope of its profile does not change signs. All 38 data profiles appear in supp. 
material Section S3. 
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(a) Gradual transition (b) Abrupt transition 




T-y x x+y t 



Figure 2: The bent-cable function, (a) A gradual quadratic transition joining two linear segments (incoming 
and outgoing). The transition period ranges from r — 7 to r + 7. Any sign change in the slope takes place at 
the critical time point (CTP). (b) An abrupt transition with 7 = yields a broken stick. The change in slope 
takes place at t, which is also the CTP. 
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10 20 30 40 50 60 

Time (in minutes) 



Figure 3: Fitted population curves (solid) with 95% pointwise credible intervals (dotted curves) for the two 
populations (grey for A and black for G). The model fit is produced assuming conditional within-individual 
independence. The estimated transition for G (i.e. M^±M y ) is marked by solid black vertical lines, and that 
for A (i.e. M TA ) by the grey vertical line. The estimated CTP for G is indicated by the dotted vertical line. 
(See supp. material Figurefflin Section S3 for the same figure but overlaid with all 38 profiles.) 
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10 20 30 40 50 60 

Time (in minutes) 

Figure 4: Observed profiles and fitted population curves (solid, in bold) with 95% pointwise credible intervals 
(dotted curves) for the two populations (grey for A and black for G). The model fit is produced assuming 
conditional within-individual independence. The estimated transition for G (i.e. m t ± M-A is marked by 
solid black vertical lines, and that for A (i.e. M TA ) by the grey vertical line. The estimated CTP for G is 
indicated by the dotted vertical line. 
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